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Using quantum Monte Carlo simulations along with higher-order spin-wave theory, bond-operator 
and strong-coupling expansions, we analyse the dynamical spin structure factor of the spin-half 
Heisenberg model on the square-lattice bilayer. We identify distinct contributions from the low- 
energy Goldstone modes in the magnetically ordered phase and the gapped triplon modes in the 
quantum disordered phase. In the antisymmetric (with respect to layer inversion) channel, the dy¬ 
namical spin structure factor exhibits a continuous evolution of spectral features across the quantum 
phase transition, connecting the two types of modes. Instead, in the symmetric channel we find a 
depletion of the spectral weight when moving from the ordered to the disordered phase. While the 
dynamical spin structure factor does not exhibit a well-defined distinct contribution from the ampli¬ 
tude (or Higgs) mode in the ordered phase, we identify an only marginally-damped amplitude mode 
in the dynamical singlet structure factor, obtained from interlayer bond correlations, in the vicinity 
of the quantum critical point. These findings provide quantitative information in direct relation to 
possible neutron or light scattering experiments in a fundamental two-dimensional quantum-critical 
spin system. 


I. INTRODUCTION 

Advances in both energy and angular resolution of 
scattering experiments enable refined experimental stud¬ 
ies of collective excitations in strongly correlated quan¬ 
tum many-body systems. A prominent example is the 
three-dimensional quantum spin-dimer system TlCuCb 
where spin excitations have been mapp ed out in great 
detail using inelastic neutron scatterin^^H. This com¬ 
pound exhibits a pressure-tuned zero-temperature tran¬ 
sition from a gapped quantum disordered phase to an 
antiferromagnetically ordered phase. In addition to iden¬ 
tifying the low-energy (transverse) Goldstone modes that 
accompany the spontaneous breaking of spin-rotation 
symmetry in the ordered phase, neutron scattering also 
detected the gapped (longitudinal) amplitude mode of 
the order-parameter field, frequently referred to as Higgs 
mode. In the three-dimensional compound TlCuGb, a 
rather successful quantitative theoretical account of these 
various spin excitation modes and the experimentally de¬ 
termined spectral weight can be obtained within a bond- 
operator mean-field description. In fact, the critical the¬ 
ory describing the underlying quantum critical point is 
the classical 0(3) field theory in four dimensions, which 
exhibits only logarithmic corrections to a Gaussian fixed 
point. This fact also implies that for nearly critical 
three-dimensional collinear antiferromagnets, the ampli¬ 
tude mode becomes increasingly sharp upon approaching 
the quantum critical poiniP. 

Upon lowering the dimensionality of the quantum mag¬ 
net, however, the effects of both thermal and quantum 
fluctuations get significantly enhanced. In particular, for 


two-dimensional Heisenberg systems, true long-range or¬ 
der is restricted to zero temperature, and the quantum 
phase transition of spin-dimer models is controlled by the 
classical Wilson-Fisher fixed point in three dimensions 
where interactions are relevaniP. These interactions in¬ 
fluence the visibility of the amplitude mode, as it can 
efficiently decay into pairs of Goldstone bosons. This 
has been examined intensively in recent years, with a 
focus towards U(l)-symmetric systems such as supercon¬ 
ductors or ultra-cold atom gases in optical lattice^SHII], 
It has been concluded from analyzing both microscopic 
models and order-parameter field theories that the am¬ 
plitude mode will be strongly masked, due to damping, 
in the order-parameter correlation function. However, 
a distinct signal of the amplitude mode can be isolated 
from accessing the so-called scalar susceptibility, in terms 
of correlations of the squared order-parameter field. For 
a magnetic systems, such as the one under considera¬ 
tion here, this implies that the amplitude mode, while 
not directly visible in the dynamical spin structure fac¬ 
tor, should be accessible via appropriate scalar corre- 
lati on f unctions, e.g., of bond-based spin-singlet opera- 
torPni. The feasibility of using probes coupling to sin¬ 
glet observables in order to detect the amplitude mode 
in (three-dimensional) magnets, e.g. via light scattering, 
has recently bee n con sidered also within bond-operator 
mean-field theorj^^E^. From field-theoretical arguments 
it is expected that a marginally damped amplitude mode 
may be detected in close vicinity of the quantum critical 
point, while deep in the ordered regime it merges with 
the continuum of multi-magnon excitations, which ex¬ 
hibits a diverging spectrum (e.g. scaling as ~ l/w with 
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the frequency uj at the ordering wave vectoi^^), mask¬ 
ing the amplitude mode. While such general conclusions 
may be drawn based on field-theory considerations, it is 
of genuine interest - also for future experimental probes 
- to provide more quantitative microscopic details on the 
excitation spectra of two-dimensional quantum magnets 
across quantum critical points. 

In this paper we use quantum Monte Carlo simulations 
combined with a stochastic analytic continuation ap¬ 
proach to access the dynamical spin structure factor and 
adequate scalar response functions for the spin-half quan¬ 
tum Heisenberg model on the two-dimensional square- 
lattice bilayer. This model has been establishecpS*^ as 
a basic spin model exhibiting a quantum critical point: 
As a function of the ratio g = J' / J between the inter¬ 
layer (J') to the intralayer (J) exchange coupling, it ex¬ 
hibits a quantum phase transition at a critical raticP^ 
Qc = 2.5220(1) that separates the small-g collinear anti- 
ferromagnet from the large-g quantum disordered dimer 
paramagnet. This model has been analysed rather inten¬ 
sively in the past, and in particular large-scale quantum 
Monte Carlo simulations have identified both the above 
quoted location of the quantum critical point and veri¬ 
fied the three-dimensional 0(3) universality class of the 
quantum phase transition. In fact, due to the absence of 
geometric frustration, this system can be studied using 
quantum Monte Carlo without any sign-problem, even in 
close vicinity of the quantum critical point, by using by¬ 
now standard cluster update algorithms. More recently, 
using an extended ensemble approach, also the quantum 
entanglement properties across the quantum phase tran¬ 
sition have been analysed and possible universal contri¬ 
butions in the bipartite entanglement have been identi- 

fiecP^. 

However, no detailed account on dynamical properties 
based on quantum Monte Carlo simulations has been 
provided thus far. Clearly, detecting the characteris¬ 
tic excitation modes of both phases in realistic spectral 
probes is of most interest here. To access such dynamical 
properties based on quantum Monte Carlo simulations, 
one requires high-quality statistical data in order to per¬ 
form the analytic continuation from the imaginary-time 
quantum Monte Carlo data to real frequencies. We note 
that dispersion relations of the excitation modes have 
been studied previously, based on various analytical ap¬ 
proaches such as bo nd-operator mean-field theory or se¬ 
ries expansion^2lIlll However, a more direct relation to 
experimental probes requires also a quantitative evalua¬ 
tion of the spectral-weight distribution. In that respect, 
one of the fundamental questions is how one crosses over 
from two gapless spin-wave excitations in the limit of 
weakly coupled planes to gapped triplon excitations at 
strong interplane coupling. This issue is of experimental 
relevance as well. For instance, the material BaCuSi 2 06 
is believed to realizd^lHIZltiie spin-half bilayer Heisenberg 
model on the square lattice (interactions between the bi¬ 
layers are weak and frustrated). These results are also 
expected to be of interest in the discussion of the bilayer 


iridate Sr 3 lr 2 07 l ^^ * ^^ l 

In the following, we provide a detailed account on the 
spectral-weight distribution in the dynamical spin struc¬ 
ture factor and also address the detection of the ampli¬ 
tude mode by appropriate scalar response functions. The 
paper is organised as follows: in Sec. H, we introduce the 
model Hamiltonian and the used computational and an¬ 
alytical methods. The dynamical spin structure factor 
of the bilayer model is then analysed in Sec. HI, while 
Sec. IV concentrates on the detection of the amplitude 
mode. The results of our analysis are summarised in the 
concluding Sec. V. Several computational details are pro¬ 
vided in the appendices. 


II. MODEL AND METHODS 

The spin-half bilayer Heisenberg model on the square 
lattice is described by the Hamiltonian 

h = s, 2 + J^(S.i • S,1 + s,2 • s,2), (1) 

* (i.i) 

where spins reside on the lower (^ = 1) and upper 
(/r = 2) layer within the i-th unit cell of a square lat¬ 
tice, where the lattice constant is set to a = 1 in the 
following. Note that each unit cell contains an interlayer 
(rung) bond of the bilayer lattice. Here, J (J') denote 
the intralayer (interlayer) exchange coupling. The model 
exhibits, in addition to the internal SU(2) spin symmetry 
and the square-lattice space-group symmetry, also a layer 
inversion symmetry in the layer indices (1 and 2). We ac¬ 
count for this additional quantum number when assign¬ 
ing a third component to an originally two-dimensional 
momentum space vector, such that k = {kx,ky,kzy, 
with fcj = 0 or TT, denoting the symmetric and antisym¬ 
metric channel with respect to layer inversion, respec¬ 
tively. Correspondingly, in position space, each spin is 
also assigned a transverse position, with respect to its 
layer index, such that is a three-component position 
vector, with the third component equal to 0 (I), for /i = I 
(/r = 2), respectively. 

Of particular interest for our analysis is the dynamical 
spin structure factor, which is defined with respect to the 
Heisenberg-picture time evolution of the spin operators 
as (with Ns denoting the number of spins) 

Ssiuj^k) = dt Y, e*(-‘-‘^'('---'--»(S.^(t).S,.(0)). 

( 2 ) 

We will also refer to the kz = 0 (kz = tt) cases as the 
symmetric or even (antisymmetric or odd) sector. In the 
presence of long-range antiferromagnetic order, one may 
also distinguish the components of 5's(w,k) parallel and 
transverse with respect to the order parameter direction, 
in which case 5'5(a;,k) represents a rotational average 
probed by the quantum Monte Carlo simulations. 
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In addition to the spin-spin correlations, we also con¬ 
sider in the following interlayer singlet bond (or spin- 
exchange) terms, 

B, = S,i • S, 2 , (3) 

which define a corresponding scalar response function in 
terms of the dynamical structure factor (with N denoting 
the number of interlayer bonds, and Ng = 2N) 

5B(cc,k) = l f (4) 

where here k and the denote two-dimensional square 
lattice k-space and lattice position vectors (i.e. with a 
vanishing third component), since the singlet operators 
Bi reside at positions on the square lattice spanned 
by the interplanar (J') rung bonds. Given its scalar na¬ 
ture, we refer to £' 3 ( 0 ;, k) also as the dynamical singlet 
structure factor. 

After having introduced the model and the relevant 
observables for our study, we next give here an overview 
of the various methods that we used. Details on the cal¬ 
culations with these methods can be found in the various 
appendices, as indicated below. 

A. Quantum Monte Carlo approach 

For the quantum Monte Carlo calculations, we used 
the stochastic series expansion method with directed loop 
updatep2H32]_ simulations, we considered finite sys¬ 

tems with Ng = 2L^ lattice sites and periodic boundary 
conditions in both square lattice directions x = (1,0)''', 
and y = (0,1)'''. In order to access ground state prop¬ 
erties, the inverse temperature /3 must be chosen suffi¬ 
ciently large. This typically requires /3J > 2L. In our 
simulations, we considered mainly finite systems with 
L = 20 with Ng = 800 sites at /3J = 50, unless speci¬ 
fied otherwise. 

To access the dynamical spin structure factor using 
the quantum Monte Carlo simulations, we efficiently^ 
measure the imaginary-time displaced spin-spin correla¬ 
tion functions directly in Matsubara frequency represen¬ 
tation, which is related to 5's(a;,k) via 

w Cl - 

S's(ia;„,k)= / duj - ^ —^S's(w, k). (5) 

Jo TT oj^+ u:i 

Here, = 27rn//3 for n = 0,1,2,... are bosonic Mat¬ 
subara frequencies. One typically needs values of n up 
to 160 to access the leading l/w^ asymptotic behaviour 
of S' 5 (ia;„,k). The numerical inversion of Eq. ([^ to ob¬ 
tain £' 5 ( 0 ;, k) from the Matsubara frequency quantum 
Monte Carlo data 5's(ia;„,k) was performed using the 
stochastic analytic continuation method in the formula¬ 
tion of Ref. [34] While such numerical analytic continua¬ 
tion methods tend to broaden inherent spectral features, 
they can still resolve a low number of excitation poles 


also from a separate continuum spectral contribution, de¬ 
pending on the quality of the imaginary-time data. 

In order to efficiently access the dynamical singlet 
structure factor S'B(w,k) of the bond-bond correlations 
from the quantum Monte Carlo simulations, we mea¬ 
sured the corresponding bond-bond correlation functions 
in imaginary-time, binned over finite-width imaginary- 
time window^^. Using an appropriate kernel for the an¬ 
alytic continuation, we then relate 5'B(a;,k) directly to 
this imaginary-time binned quantum Monte Carlo data, 
without the need for an intermediate unfolding of the bin- 
resolution correlation function. Further details of this 
measurement setup are provided in App. 


B. Spin-wave theory 

Upon expressing the spin fluctuations about the clas¬ 
sical magnetically ordered state using the standard 
Holstein-Primakoff representatiorP^, one arrives in the 
harmonic approximation at the linear spin-wave descrip¬ 
tion of the magnetic excitations in the ordered phase. 
As discussed in Sec. Ill, the comparison between the 
spin-wave results and the dynamical structure factor 
as obtained from the quantum Monte Carlo simulation 
improves upon taking into account corrections beyond 
the harmonic approximation in the spin-wave expansion. 
This has been done here by considering the next-to- 
leading order in the 1/S expansion of the Hamiltonian, 
and by performing a mean-field decoupling of the bosonic 
interaction terms to renormalize the coupling constants. 
Details on both the linear and higher-order spin-wave 
theory calculations are presented in App. 


C. Perturbation theory m 1/g 

In the quantum-paramagnetic phase, an efficient de¬ 
scription of the spin dynamics can be obtained from per¬ 
turbation theory in the intralayer coupling J, i.e., start¬ 
ing from the limit of decoupled dimers which form a 
product state of singlets in the limit J = 0. Exciting 
a single rung gives rise to a triplet excitation, which can 
delocalize for finite J, giving rise to a gapped excita¬ 
tion mode. We perform a systematic expansion for the 
dynamical spin structure factor in 1 /gr up to second or¬ 
der for the triplet dispersion and to linear order in the 
spectral weight. While for the antisymmetric channel 
closed explicit expressions are obtained, we need to per¬ 
form a numerical diagonalization of the effective interac¬ 
tion Hamiltonian in the two-triplon sector to access the 
symmetric channel of the dynamical spin structure fac¬ 
tor. Details on these calculations are provided in App. 0 
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D. Bond-operator-based 1/d expansion 

Bond operator^^ were initially introduced to effi¬ 
ciently describe the paramagnetic phase of spin-half 
coupled-dimer magnets, the square lattice bilayer con¬ 
sidered here being a prominent example. In this ap¬ 
proach one introduces bosonic operators corresponding 
to spin-1 e xcitat ions (often called triplons) atop a singlet 
backgrouncP^JSH Generalisations of the bond-operator 
approach to magnetically ordered stated and to arbi¬ 
trary spirP^ have been in vestig ated in the past. More¬ 
over, it was recently showiP^ED that bond operators en¬ 
able a controlled description of coupled-dimer systems 
across the entire phase diagram using l/d as a small pa¬ 
rameter, where d is the spatial dimension: relevant ob¬ 
servables can be obtained in a systematic 1/d expansion 
once the dimer lattice has been generalized to d space 
dimensions (for details we refer to Refs. HOI and I^Tl) . In 
particular, the dispersion relation of triplon excitations in 
the disordered phase as well as that of gapless transverse 
modes and the gapped longitudinal mode in the ordered 
phase were reported along with their corresponding spec¬ 
tral weights in the dynamic spin susceptibility. As will 
be discussed in Sec. Ill, a comparison with the disper¬ 
sions obtained from quantum Monte Carlo shows very 
good agreement after a proper mapping of the model pa¬ 
rameters (the coupling ratio g), as explained in Sec. III. 
In order to study the amplitude (Higgs) mode, we also 
calculated the bond-bond correlation to leading order in 
1/d. A comparison with quantum Monte Carlo results 
will be discussed in Sec. IV. We find that to leading order 
in l/d, the longitudinal mode gives rise to a single mode 
spectral contribution to the inter-layer singlet dynamical 
structure factor. Further details on the bond-operator 
calculations can be found in App. 

III. DYNAMICAL SPIN STRUCTURE FACTOR 

The dynamical spin structure factor S's(a;,k) is di¬ 
rectly related to the scattering intensity as probed by 
inelastic neutron scattering experiments. Here, we moni¬ 
tor the evolution of Ss{uj, k) upon varying the interaction 
ratio g across the quantum critical point. 

A. Quantum Monte Carlo results 

We first present the results from the quantum Monte 
Carlo simulations for the dynamical spin structure factor. 
In Fig. the dynamical spin structure factor Ss{oJ, k) is 
shown for different values of g, each along a standard- 
path in the two-dimensional Brillouin zone, for the sym¬ 
metric {kz = 0) and the antisymmetric {kz = tt) channel, 
separately. The structure factor in the antisymmetric 
channel is dominated by a sharp single-mode-like con¬ 
tribution that softens in the antiferromagnetic regime at 
k = Q, the magnetic Bragg peak position Q = (tt, tt, tt)'''. 


while it exhibits a fully gapped branch in the quantum 
disordered regime, with a minimum gap Arig) in the dis¬ 
persion, located at k = Q. In the antiferromagnetic re¬ 
gion, the dispersive feature traces the dispersion relation 
of the antiferromagnetic magnon (Goldstone) excitation¬ 
mode, while in the quantum disordered region, it follows 
the gapped triplon dispersion. Both these relations will 
be substantiated by a direct comparison to spin-wave the¬ 
ory and perturbation theory calculations in 1/g for the 
antiferromagnetically ordered and quantum disordered 
regime respectively, as discussed below in Sec. |HIB| 

Approaching the quantum critical point, we observe 
enhanced finite-size effects due to the algebraic growth of 
the correlation length within the quantum critical region. 
We anticipate the most pronounced finite-size corrections 
to appear right at the quantum critical point: based on 
the relativistic invariance with a dynamical critical expo¬ 
nent z = 1, at low energies finite-size corrections of the 
peak position proportional to 1 /L dominate at criticality. 
Indeed, enhanced finite-size corrections at the quantum 
critical point were reported recently for the low-energy 
dispersion of the Goldstone mode^, and are visible in 
Fig.g for g = 2.522, where a small finite-size gap at 
k = Q can be clearly resolved. A finite-size study of the 
dynamical spin structure factor, shown in Fig. indeed 
reveals a linearly vanishing finite-size excitation gap as a 
function of 1/A, as shown in the inset of that figure. 

In addition to the sharp, single-mode feature, the 
structure factor also exhibits a broader distribution of 
spectral weight for energies above the single-mode thresh¬ 
old. This contribution will be discussed in more detail in 
Sec. VI, since it is of particular interest with regards 
to the possibility of observing a Higgs peak in the spin 
structure factor. 

The symmetric channel exhibits in the antiferromag¬ 
netically ordered phase a similarly sharply resolved dis¬ 
tribution of the spectral weight, which however softens 
near the F = (0,0,0)''' point. Its dispersion follows to a 
very good accuracy that of the antisymmetric mode with 
a shift of the in-plane momentum by (tt, tt). At the F 
point, the spectral weight is completely suppressed for 
all values of g: this relates to the fact that spin fluc¬ 
tuations at the F point, i.e. in the uniform magnetisa¬ 
tion, vanish, due to the SU(2) symmetry of H. The gap 
at k = (tTjTt, 0)''' decreases upon approaching the limit 
of decoupled layers (g —>■ 0), and in the limit of decou¬ 
pled layers, the structure factor in the symmetric channel 
softens at k = (tt, tt, 0)''', which is the single layer Bragg 
peak position. By contrast to the antisymmetric chan¬ 
nel, the overall spectral weight in the symmetric channel 
is strongly suppressed for large values of g: well in the 
quantum disordered regime, only a faint spectral weight 
distribution is present in the symmetric channel, while 
the antisymmetric channel provides a direct trace of the 
dispersion relation of the gapped triplon excitation mode. 

The bilayer Heisenberg model is obtained as the effec¬ 
tive low-energy theory for the spin dynamics within the 
Mott insulator, in the strong-coupling limit of the half- 
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FIG. 1. Dynamical spin structure factor Ss(a;,k) in the symmetric (upper panel) and antisymmetric (lower panel) channel for 
the spin-half Heisenberg model on the square lattice bilayer at different ratios g — J'/ J of the interlayer (J') to intralayer (J) 
coupling strength along the indicated path in the two-dimensional Brillouin zone. 



FIG. 2. (Color online) Dynamical spin structure factors 
^^(aijk = Q) for the spin-half Heisenberg model on the 
square lattice bilayer at the quantum critical point, g = g^. 
The finite-size scaling of the excitation gap is shown vs. 1/L 
in the inset. 


filled Hubbard model on the square lattice bilayer. At 
strong local Hubbard repulsion U t, we thus expect 
that the evolution of the spectral weight in the dynami¬ 
cal spin structure factor, as described above, can be ob¬ 
served also in the Hubbard model, upon tuning the ratio 
t' jt of the interlayer tunneling (t') to the intralayer tun¬ 
neling it). It is less clear, however, that this also holds 
in the intermediate coupling regime, where U is of order 
the bandwidth, and where residual charge fluctuations al¬ 
low for higher-order spin exchange processes. To address 
this question, we performed also determinantal quantum 


Monte Carlo simulations for the Hubbard model on the 
square lattice bilayer. From our simulations, we find that 
also in the intermediate coupling region the dynamical 
spin structure factor exhibits the characteristic behaviour 
that we observed in the Heisenberg limit. We provide de¬ 
tails of these results for the Hubbard model in the supple¬ 
mental material to this papei!^. Returning here to the 
Heisenberg model, we next perform a detailed compari¬ 
son of the quantum Monte Carlo results to the various 
theoretical approaches that we listed in Sec. H. 


B. Comparison to analytic results 

In the following, we compare our quantum Monte Carlo 
data to the results obtained within (i) linear and higher 
order spin-wave theory, (ii) the l/g perturbation the¬ 
ory as well as (iii) the 1/d expansion within the bond- 
operator approach. While the spin-wave and the large-(/ 
approach are restricted to the ordered and quantum dis¬ 
ordered regimes, respectively, the 1/d expansion allows 
us to calculate appropriate structure factors throughout 
the whole phase diagram. The upper two rows of Fig. 
and the upper row of Fig. show the dispersion rela¬ 
tions for the region below and above the quantum crit¬ 
ical point, respectively. In these panels, the mode dis¬ 
persions from the quantum Monte Carlo data have been 
obtained by the peak positions in the spectral weight dis¬ 
tribution, with the symbol size indicating the estimated 
uncertainty. In addition to comparing the dispersions, 
we also compare, in the lower rows of Fig. [^and Fig. 
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FIG. 3. (Color online) Comparison of the dynamical spin 
structure factor at different ratios g = J'/J inside the an- 
tiferromagnetically ordered region for the spin-half Heisen¬ 
berg model on the square lattice bilayer between the quantum 
Monte Carlo (QMC) results and linear and 2nd order spin- 
wave theory (SWT) as well as bond-operator theory (BOT), 
both in harmonic approximation and including leading 1/d 
corrections. The upper panel compares the dispersion of the 
single-mode contribution, and the lower panel the integrated 
spectral weight ^^(k). Here, the results for the symmetric 
(antisymmetric) channel have been normalised by the value 
of Ss{{Tf, TT, 0)) ( 53 ( 0 )), corresponding to the maximum (min¬ 
imum) value in that channel. 

the integrated spectral weight 

SsCk) = J duj SsiuJ,^), (6) 

which in the quantum Monte Carlo simulations is con¬ 
veniently obtained from the equal-time spin-spin corre¬ 
lations, to the results from spin-wave theory and the 1 /g 
perturbation theory as well as the bond-operator the¬ 
ory. For the ordered phase, we performed the bond- 
operator theory calculations for the integrated spectral 
weight only in the harmonic approximation, while for the 
disordered phase, we also considered the leading 1/d cor¬ 
rections. We also note that, in contrast to the mode dis¬ 
persion, the integrated spectral weight can be obtained 




{kx,ky,ki) {kx,ky,kx) 

FIG. 4. (Color online) Comparison of the dynamical spin 
structure factor in the antisymmetric channel at different ra¬ 
tios g = J'/J well inside the quantum disordered region for 
the spin-half Heisenberg model on the square lattice bilayer 
between the quantum Monte Carlo (QMC) results and 1/g 
perturbation theory (PT) as well as bond-operator theory 
(BOT) in both harmonic approximation and with leading 1/d 
corrections included. The upper panel compares the disper¬ 
sion of the single triplon contribution, and the lower panel the 
integrated spectral weight ^^(k). 


directly from the quantum Monte Carlo simulations with¬ 
out the need of performing the analytic continuation. 

We first compare our numerical results to the disper¬ 
sions obtained from spin-wave theory. From the upper 
two panels of Fig. we find that (i) linear spin-wave 
theory, while providing a good overall account of the spin- 
wave dispersion, systematically underestimates the spin- 
wave energies, and (ii) the higher order spin-wave theory 
approximation provides a significant improvement to the 
overall dispersions. As detailed in App. the net effect 
of the second order corrections to linear spin-wave theory 
can be expressed in terms of a (/-dependent renormaliza¬ 
tion of the coupling ratio g, which for g > 1 leads to an 
enhanced effective coupling ratio geff > g, which results 
in the hardening of the spin-wave dispersion, as observed 
in Fig. As discussed in App. we find that an even 
further, heuristic renormalization of g^s allows us to ob¬ 
tain even better matches of the effective spin-wave theory 
dispersion to the numerical results. 

With respect to the integrated spectral weight 5's'(k), 
which is shown in the lower two panels of Fig. we 
find that up to intermediate coupling ratios (cf. the case 
g = 1 in the left panels of Fig. |^, both linear and second 
order spin-wave theory provide a good overall account of 
the spectral weight distribution in the symmetric chan¬ 
nel. Upon closer inspection however, one notices that 
the spin-wave theory results for g = 1 fall slightly be¬ 
low the quantum Monte Carlo data for S's(k) for k near 
(tt, 0,0)'''. This difference gets more pronounced upon in¬ 
creasing g towards the quantum critical point, cf. the 
data for g = 2 in the right panels of Fig.|^ The failure of 
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low-order spin-wave theory to fully account for the spec¬ 
tral properties for wave vectors near k = (tt, O)’’’ has been 
noticed previously for the case of a single layer Heisen¬ 
berg model (i.e., g = O', cf., e.g., Ref. HHfor an extended 
discussion) and was recently linked to strong attractive 
magnon-magnon interactions that induce a Higgs reso¬ 
nance from two-magnon state^^. One would then expect 
that the increase in the deviation to low-order spin-wave 
theory that we observe in the integrated spectral weight 
for g = 2 relates to an enhancement in the Higgs res¬ 
onance formation. If fact, as discussed in Sec. HI, we 
can identify a well defined low-energy Higgs peak in the 
singlet structure factor for g near gc, which adds further 
support to this scenario. 

For the antisymmetric mode, we find that both linear 
and higher-order spin-wave theory accounts well for the 
integrated spectral weight distribution for g = 1, while 
they overestimate the integrated spectral weight near the 
ordering wave vector, k = Q, closer to gc, cf. the data 
for g = 2 in Fig. reflecting the fact that quantum 
fluctuations of the order parameter are underestimated 
within these approximations. 

Next, we focus on the quantum disordered region, and 
discuss the comparison to the 1/g perturbation theory. 
Well inside the quantum disordered region, we find that 
the 1/g perturbation theory, presented in App. pro¬ 
vides a rather good account of the quantum Monte Carlo 
data. A corresponding comparison of the triplon disper¬ 
sion and the corresponding spectral weight (i.e. in the 
antisymmetric channel) for g = 6 and g = 10 is shown 
in Fig. While the triplon dispersion is obtained very 
accurately within our 1/g perturbation calculations for 
both values of g, the integrated spectral weight at g = 6 
already exhibits somewhat larger quantitative differences 
than at g = 10. As detailed in App.[C| we performed the 
perturbative calculations for the triplon dispersion up to 
quadratic order in 1/g, while the spectral weight was 
obtained up to linear order, and thus is expected to be 
less accurate at larger values of 1/g. Nevertheless, the 
overall shape of the spectral weight distribution is well 
represented at this leading order already. 

Thus far, in discussing the quantum disordered region, 
we considered only the antisymmetric channel, which 
is dominated by the gapped triplon mode with a pro¬ 
nounced spectral weight. However, within 1/g perturba¬ 
tion theory, we can also obtain a quantitative description 
of the faint spectral weight distribution for the symmet¬ 
ric channel in the quantum disordered regime. Figure 
shows the result for the dynamical spin structure factor in 
the symmetric channel as obtained from the 1/g pertur¬ 
bation expansion by a numerical diagonalization of the 
effective Hamiltonian within the two-triplon sector for a 
L = 10 system. The structure factor is seen to display 
a concentration of the spectral weight at k = (tt, tt, 0)''', 
which corresponds to an eigenstate of the effective Hamil¬ 
tonian in that sector. This is in good agreement with the 
quantum Monte Carlo data that also exhibit a concen¬ 
tration of the spectral weight for the symmetric chan- 



{/k^, ky, k^) 

FIG. 5. (Color online) Distribution of the amplitude in the 
dynamical spin structure factor 53 ( 0 ;, k) for the symmetric 
sector (kz = 0) as obtained within the 1/g perturbation the¬ 
ory for the numerical diagonalisation of the effective Hamil¬ 
tonian within the two-triplon sector on a L = 10 system. The 
sizes of circles are proportional to the amplitude at the centre 
positions. Dashed lines indicate the upper and lower bounds 
of the continuum. 

nel near k = (tt,tt, O)’’’, cf., e.g., the data for g = 3 in 
Fig.[T} One furthermore finds from the perturbative cal¬ 
culations, that the amplitude of the symmetric channel is 
of order (1/g)^, and thus vanishes in the limit J/ J' —)■ 0 
(cf. App. for details of the calculation), which is re¬ 
flected also by the quantum Monte Carlo data. 

The comparison of the quantum Monte Carlo data to 
the bond operator-based 1/d expansion requires an ap¬ 
propriate mapping from the coupling ratio g to the pa¬ 
rameter q (defined as q = Jd/ J' = d/g, with d the dimen¬ 
sion of the system, here, d = 2) that enters the 1/d ex¬ 
pansion. Indeed, within the 1/d expansion, the quantum 
critical point is obtained as qc = l/2-|-3/(16d)-|-Cl(l/d^), 
so that at the harmonic level (i.e. without 1/d correc¬ 
tions) qc = 1/2, while including first-order corrections in 
1/d leads to qc = 0.59375, to be compared to the value of 
{gQMc)c = d/gc = 0.7930(2) that results from the (accu¬ 
rate) quantum Monte Carlo estimate of gc- We need to 
account for this difference when performing the compar¬ 
ison, in particular in the vicinity of g^ For a given value 
oi g < gc, we fix the value of q such that {q — qc) equals 
the corresponding absolute distance in qqMC = d/g to 
{qQMc)c, i-e., such that 

g- qc = qQMc - iqQMc)c- (7) 

In the quantum disordered phase, g > gc, it is more ap¬ 
propriate to relate the parameters by considering a hxed 
relative distance to the quantum critical point, i.e., a 
value of q is obtained, such that 

q- qc _ qQMc - {qQMc)c 

qc {qQMc)c 







The corresponding mode dispersions are shown in the up¬ 
per panels of Fig. |^for g < gc and Fig. 0 for g > gc, re¬ 
spectively, while the lower panels of these figures show the 
comparison in the integrated spectral weight <S's(k). We 
find that with the above parameter mappings, the bond- 
operator theory provides a good qualitative account of 
both the spin-wave dispersion and the triplon dispersion 
within this unified approach, at least when the leading 
1/d corrections beyond the harmonic approximation are 
taken into account. Nevertheless, the harmonic approxi¬ 
mation already keeps track of the leading dispersive fea¬ 
tures and provides a good account of the integrated spec¬ 
tral weight in both phases. For the following discussion of 
the amplitude mode, we thus limited the bond-operator 
theory calculations to the harmonic level. 

IV. AMPLITUDE (HIGGS) MODE 

The longitudinal, amplitude fluctuations of the order 
parameter field need to be properly taken into consider¬ 
ation in a quantitative description of the quantum crit¬ 
ical bilayer He isenb erg antiferromagnet and its low en¬ 
ergy propertie^^nmi Por example, within bond-operator 
theory calculations, the amplitude mode is found to 
soften upon approaching the quantum critical point, thus 
restoring the SU{2) symmetry, while it hardens towards 
an elevated energy scale of order 4J in the limit of de¬ 
coupled layers {J' = 0), rendering it insignificant with 
respect to the or der-p arameter fluctuations well inside 
the ordered phase^SEll While the relevance of the am¬ 
plitude mode on e.g. the order parameter strength and 
the critical value of g has been well studied in the past, 
only recently have possible direct signatures of the ampli¬ 
tude mode and appropriate experimental probes been ad¬ 
dressed. In particular, it has been argued that response 
functions which couple to the square of the order param¬ 
eter will show a distinct peak from the amplitude mode 
(Higgs peak), well separated from the multi-magnon con¬ 
tribution to the overall spectral weight in t he ordered 
phase, close to the quantum critical pointPi^. 


A. Quantum Monte Carlo results 

To quantify these considerations for the bilayer Heisen¬ 
berg model, we examine in this section the interlayer- 
bond spin-exchange correlations, in terms of the corre¬ 
sponding dynamical singlet structure factor S'B(w,k) in¬ 
troduced in Sec. H. Based on its scalar nature, we expect 
the correlations of this observable to exhibit in addition 
to possible multi-magnon contributions also a distinct 
signal from the amplitude mode, most pronounced in the 
low-energy regime, and within the vicinity of the quan¬ 
tum critical point. In Fig. we show our numerical 
data for different values of g, taken on a L = 20 system 
at PJ = 50, representative of ground state expectation 
values. We obtain for all values of g a spectral contri¬ 


bution to S'b(w, k) at elevated energies uj « 4J, which 
does not exhibit a strong dependence on the coupling 
ratio g. In addition however, we observe in the vicin¬ 
ity of the quantum critical point a distinct feature in 
the low-energy region, i.e., separated from the higher- 
energy signal. This additional spectral branch is most 
pronounced in the vicinity of the F point, where also its 
spectral weight is most clearly visible. 

To better exhibit the emergence of the low-energy 
Higgs peak from the multi-magnon contribution, we show 
in Fig.j^the dynamical singlet structure factor Sb{u}, k = 
0) at the F point for different values of the coupling ra¬ 
tio gc- Within the range “2 < g < gc, detect a sep¬ 
arate, low-energy peak split-off from the broad second 
peak in the spectral weight at elevated energies. The 
position of the low-energy peak tends towards lower en¬ 
ergies upon tuning g towards its critical value gc- How¬ 
ever, approaching the quantum critical point, we again 
observe enhanced finite-size effects, most pronounced at 
the quantum critical point itself. This behavior is illus¬ 
trated for g = gc va Fig. We find that for g below 
about g « 2, the Higgs peak merges with the second 
peak, so that a Higgs peak cannot be discerned anymore 
as a separate spectral feature, restricting its observation 
to the close vicinity of the quantum critical point. 

B. Comparison to bond-operator theory 

Within the parameter regime for which we can identify 
a Higgs-mode peak at the F point, we also observe traces 
of the corresponding mode dispersion near F. It is thus 
interesting to compare our quantum Monte Carlo data to 
the amplitude-mode dispersion obtained from the bond- 
operator-based I /d-expansiorP3. Referring for details on 
the 1/d-expansion calculations oi Sb to App.|^ we com¬ 
pare in Fig. the resulting dispersion relation for the 
amplitude mode with the quantum Monte Carlo data, 
using the parameter mapping from Eq. ([^. 

Using bond-operator theory we also calculated the 
spectral weight of the amplitude mode within the har¬ 
monic approximation. As discussed in more detail 
in App. one finds an enhanced spectral weight in 
S'B(a;,k) near the F point, in particular in the vicinity 
of the quantum critical point, a result which appears to 
agree with the quantum Monte Carlo data, which also 
exhibit more pronounced spectral weight for the ampli¬ 
tude mode branch near the F point and near the quantum 
critical point. A more quantitative comparison is limited 
by the fact that there are contributions (of multi-mode 
character) to the spectral weight observed in beyond 
the pure Higgs-mode contribution considered in App. 

C. Scaling form and comparison to S's(aj,k) 

To analyse further the shape of the Higgs peak, we 
compare our numerical results to the universal low- 
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FIG. 6. (Color online) Dynamical singlet structure factor 5's(a;,k) for the spin-half Heisenberg model on the square lattice 
bilayer at different ratios g = J'/J of the interlayer (J') to intralayer (J) coupling strength along the indicated path in the 
two-dimensional Brillouin zone. Also indicated by lines are the amplitude-mode dispersions obtain within the bond-operator 
theory. 





FIG. 7. (Color online) Dynamical singlet structure factor 5s(u;,k = 0) at the F point for the spin-half Heisenberg model on 
the square lattice bilayer at different ratios g = J'/J of the interlayer (J') to intralayer (J) coupling strength. 


energy scaling form of the scalar response functiorP ob¬ 
tained from a 1/A^-expansion of the 0(N) model in the 
quantum critical region, 

S{uj) oc A3-2 /-$(w/A), (9) 

in terms of the characteristic energy scale in the quan¬ 
tum critical region, A oc — g/JcY^ with the correlation 
length exponent v of the three-dimensional 0(3) univer¬ 
sality class, to which the quantum critical point gc in the 
bilayer Heisenberg model belongs. The value of the crit¬ 
ical exponent v has been determined by previous large 
scale Monte Carlo simulations^ as v = 0.710(2). Fi¬ 
nally, $ denotes a corresponding scaling function. We 
estimate the (/-dependent excitation gap of the ampli¬ 
tude mode ruH (g) (Higgs mass) from the position of the 
low-energy peak in iS'B(w,k = 0 ), and the energy scale 
A{g) for g < gc in terms of the triplet excitation gap 
Arig') in the quantum disordered phase, where g' > gc is 
the mirrored (with respect to the quantum critical point) 
coupling ratio defined by g' — gc = Jc — g- In other words. 


A{g) = Arigc + {gc - g))- For 2 < 5 < gc, i.e. in the 
region where we can identify a Higgs peak, we find that 
mH{g)/A{g) = 2.6(4), in agreement with a previous es¬ 
timate for this ratio from the effective 0(3) field theorjP. 

Based on this relation, we show in Fig. I^the appropri¬ 
ately rescaled data for Sb{oj, k = 0 ) from different values 
of g according to the scaling law in Eq. (©• From this 
analysis we find that the low-energy part of iS'_b(w, k = 0) 
relates well to this scaling form, adding further support 
to the interpretation of the low-energy peak in terms of 
the amplitude mode. By contrast, the second peak at en¬ 
ergies a; « 4 J does not follow this scaling form, hence is 
not part of the universal signal in the response function. 
In fact, as already mentioned above, we see from Fig. 
that the position of the second peak does not change 
much upon varying g, in contrast to the Higgs peak. 
The observation, that the amplitude mode follows the 
scalar response function scaling form also implies that in 
this two-dimensional system, the amplitude mode is only 
marginally damped in the vicinity of the quantum crit- 
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FIG. 8. (Color online) Dynamical singlet structure factor 
= Q) for the spin-half Heisenberg model on the 
square lattice bilayer at the quantum critical point, g — Qc- 
The finite-size scaling of the excitation gap is shown vs. 1/L 
in the inset. 



FIG. 9. (Color online) Scaling plot of the dynamical singlet 
structure factor ^^(a;, k = 0) at the F point for the spin-half 
Heisenberg model on the square lattice bilayer at different 
ratios g = J'/ J of the interlayer (J') to intralayer (J) coupling 
strength near the quantum critical pint g^. 


ical point, i.e. with a width that scales like the Higgs 
mass niH near the quantum critical point, thus making 
the observation of this excitation mode possible. 

We find that the quantum Monte Carlo data are consis¬ 
tent also with the low-energy scaling predictioiP S'(w) oc 
in as much as the imaginary-time (t) data for the 
bond-bond correlations at k = 0 (which are actually cal¬ 
culated in the simulation) agrees with the corresponding 
decay Sb (t, k = 0) oc 1 /r^ at large imaginary times (cf. 
App. a for details). However, we cannot independently 
estimate the large-r scaling behaviour since we are lim¬ 
ited by enhanced statistical noise on the imaginary-time 
data in the relevant r regime (cf. App. 0 . 

In order to contrast further the appearance of the am¬ 
plitude mode in the dynamical singlet structure factor 


as compared to the dynamical spin structure factor, we 
show for different values of g in Fig. 10 both (i) the dy¬ 
namical singlet structure factor 5'b(u;, k = 0 ) at the wave 
vector where Sb softens at gc, and (ii) the dynamical spin 
structure factor Ss{uj, k = Q) at the Bragg peak position 
where a possible Higgs mode-contribution would soften 
at the quantum critical poiniPil. 

As mentioned already in Sec. |mg in the ordered 
phase, Ss{uj,]i) exhibits in addition to the magnetic 
Bragg peak an extended tail of spectral weight. This 


is shown for k = Q for several values of g in Fig. 10 


Atop a background signal that falls off with increasing 
w, and which is expected from general considerations in 
the symmetry-broken phas^^^, we observe two broad ad¬ 
ditional features in Ss{uj,Q) at w « J and ui « 4J, re¬ 
spectively. Compared to the magnetic Bragg peak, both 
these features are significantly broader and of suppressed 
spectral weight. These features do not exhibit system¬ 
atic trends upon varying g, apart from weak shape vari¬ 
ations for different values of g. By contrast, the data 
for £' 3 ( 0 ;, k = 0) exhibits a clear Higgs peak split-off 
at g « 2 , that we already analysed above, with a g- 
dependent peak position and a shape that agrees with the 
scalar response function scaling form. Fig. [T 0 | illustrates 
how the presence of the extended residual spectral weight 
tail beyond the magnetic Bragg peak masks a possible ob¬ 
servation of the low-energy Higgs peak in the dynamical 
spin structure factor even near criticality, in contrast to 
the scalar response function, which is strongly suppressed 
at low energies. Furthermore, we find that both structure 
factors show an enhanced spectral weight in the form of 
a broad peak at similarly elevated energies w « 4 J, sug¬ 
gesting that the spectral weight observed in ^^(Wjk) in 
this energy range is dominated by multi-magnon contri¬ 
butions. 

We finally note that we also examined using the 
quantum Monte Carlo simulations the intraplanar bond- 
bond correlations, and find in the symmetric sector 
(with respect to layer inversion) a similar (but weaker) 
amplitude-mode contribution emerging in the low-energy 
region, while no such feature is observed for the antisym¬ 
metric channel. As noted in App. this observation 
also agrees with the bond-operator-based 1 /d expansion 
calculations of the intraplanar bond-bond correlations. 


V. CONCLUSIONS 

Using quantum Monte Carlo simulations and a vari¬ 
ety of analytical approaches, we have obtained precise 
information on the dynamical response functions of the 
spin-1/2 Heisenberg model on a bilayer system of square 
lattices. We have shown that the transition from two 
gapless spin-wave excitations in the limit of weakly cou¬ 
pled planes to a single gapped triplon excitation in the 
strong coupling limit occurs via a rapid suppression of the 
spectral weight of the symmetric spin-wave mode upon 
increasing the interlayer coupling. The intensity of that 
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FIG. 10. (Color online) Comparison between the dynamical spin structure factor ~ Q) and the dynamical singlet 

structure factor 5s(a;,k = 0) for the spin-half Heisenberg model on the square lattice bilayer at different ratios g = J'/J of 
the interlayer (J') to intralayer (J) coupling strength. 


mode in the dynamical spin structure factor probed in 
neutron scattering is already extremely weak at the quan¬ 
tum critical point where the system becomes gapped, 
and it is the antisymmetric mode that smoothly con¬ 
nects to the triplon mode of the strong coupling limit. 
We have also looked for traces of the amplitude (Higgs) 
modes in various dynamical correlation functions. This 
mode is clearly visible as a distinct spectral feature in the 
interlayer-bond dynamical singlet structure factor, which 
opens to the way to its detection with light scattering. 
However, this mode is not visible in the dynamical spin 
structure factor, most likely because it is masked by an al¬ 
most dispersionless continuum of incoherent excitations. 
It is our hope that these conclusions will further motivate 
the experimental investigation of the dynamical response 
of bilayer systems. 
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Appendix A: Bond operator measurements in 
quantum Monte Carlo 


and B that are constructed from Hamiltonian terms, such 
as the Fourier modes 

= (Al) 

of the spin-exchange terms of the bond-operators Bi in¬ 
troduced in Sec. H, instead of the imaginary-time (r) 
correlation function {A(t)B), we measure an auxiliary 
correlation function that is defined in terms of a discrete 
time grid laid over the t- interval [0,/3], with width St. 
Namely, we consider the Monte Carlo estimator 

Cab{1 5t) = NAe[iST,(i-i)ST] NBe[ST,o]^ , (A2) 

where (•)mc denotes the Monte Carlo expectation value, 
and Na^i counts the number of operators A in the SSE 
operator sequence which lie within the time imaginary- 
time interval I, while I = with L^. given by the 

number of time-slices, /3 = Lt St. Using the periodicity 
of the SSE configuration, one can further improve the 
above estimator to 

CAsil 5 t) = — {N(^a,b),i)mc ’ (A3) 

where N(^a,b),i counts the number of pairs {A, B) in the 
operator string where A acts I — 1 time intervals later 
than B. This auxiliary correlation function is related to 
the correlation function {A(t)B) through a convolution, 

Cab{ISt) = J Q{St-\t\){A{{1-1)St+t)B). 

.(A4) 

Instead of performing an approximate de-convolution of 
this relatioiP^, here we make use of the convoluted kernel 


To efficiently calculate the dynamical singlet struc¬ 
ture factor S'B(a;,k) within the stochastic series expan¬ 
sion (SSE) formulation of quantum Monte Carlo simu¬ 
lations, we extend an approach put forward in Ref. 1331 
It is based on the well-known mapping of the discrete 
SSE configuration space onto a continuous-time world¬ 
line formulatioiP^. For hermitian conjugate operators A 


Ksr{uJ,T) 



which yields 


(AS) 


Ksr{u},T) 


(giro; _ e(2r-/3)u;) 


(A6) 
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S[ = , 

Sf = Sf , (B2) 

Sz = e^Q '^'Sf . 

In this rotated frame, the classical ground state is thus 
ferromagnetic by construction. 

0 1 2 

t/J ^ 1. Linear spin-wave theory 



FIG. 11. (Color online) Imaginary-time (r) dependence of 
S's('r, k) = (i3k('7")B_k) at the F point, k = 0, for a cou¬ 
pling ratio of p = 2.1, along with an asymptotic l/r‘*-decay, 
indicated by the line. 


in order to directly relate the quantum Monte Carlo es¬ 
timator Cab{1 St) to Sab{oj) = f dt corre¬ 

sponding to the real-time correlation function {A{t)B), 

poo 1 

Cab{15t)= J -^Ksr{uj,T)SAB{to). (A7) 

This relation holds true without any approximation, ir¬ 
respectively of the discretisation used in setting up the 
Monte Carlo estimator. The final inversion problem is 
then solved using e.g. the stochastic analytic continua¬ 
tion approach, with the kernel for the present 

case. As an example, we show in Fix. the quantum 
Monte Carlo data for k) = {B]^{t)B^\^ at the T 
point, k = 0, and a coupling ratio oi g = 2.1, along with 
the indicated asymptotic 1/r'^-decay. 


Appendix B: Spin-wave theory 

In this appendix, we provide details on both the lin¬ 
ear spin-wave theory calculations that we performed on 
the square-lattice bilayer as well as on the higher-order 
spin-wave theory approach that we used. In the classical 
limit of large spin S', the ground state of the Heisenberg 
model on the bipartite square lattice bilayer consists of 
an antiferromagnetic arrangement of the spins within and 
between the layers, irrespectively of the coupling ratio 
g. The classical spin structure can be parametrized by 
means of a single helix with pitch vector Q = ( 7 r, 7 r, 7 r)''', 
such that for a spin on site I, 

Si = S(0,0, . (HI) 

In this appendix, the index I runs over all sites of the 
lattice, where I = {i,g) is given by a rung index i plus 
a layer index y, as in the main text. The quantum spin 
Hamiltonian is now rewritten expressing the spin opera¬ 
tors in the local basis of the classical spin orientations, 


In the quantum model, the deviations from the classi¬ 
cal order are expressed via the Holstein-Primakoff repre¬ 
sentation of spin operator^^. To first order in 1/S, the 
expressions take the form 


Sf = 


Sf 


Sf = 




(o; -f aj) 


I ) ’ 


V25 t, 
— (o; - a)), 


(B3) 


2i 
S — n 


i ’ 


with (a|) local bosonic annihilation (creation) oper¬ 
ators, and ni = a|a;. In this bosonic representation 
the expression of the quantum Hamiltonian, truncated 
at the harmonic order, takes the following compact form 
in Fourier space 



where Ns is the total number of lattice sites, and the 
coefficients Ak and i?k are given by 


j' 

Ak = 2J -|- 

J' 

i?k = — J(cOsfca; -I- cosky) - — coskz, 


(B5) 


with kz = 0 {kz = tt) for the sym metr ic (antisymmetric) 
sector. Terms of order S^ in Eq. (B4) correspond to the 


classical energy and terms of ord er S to the 1/S correc¬ 
tion. The quadratic Hamiltonian (B41 is diagonalized by 


introducing new bosonic quasiparticle operators via 
“k = + ^k«U> 

4 = ^ik4 + ^k«-k> 


(B6) 


with Uk and Wk given by 


Uk = 


-Sgn(Hk) 


I 




+ 1 


Vk = 




- 1 


(B7) 
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The diagonal representation of H reads 

H = -SiS + 1) (^2J|| + iV5 + ^ > 

(B8) 

with the free magnon dispersion given by 

Wk = 2S^Al-Bl (B9) 


The ground state of (B8| corresponds to the vacuum |0) 
of the Bogolyubov quasiparticles a and, for S = 1/2, 
we obtain a symmetric and an antisymmetric mode with 
dispersions 


= 2T^[l - -i{k,;,ky)\ [l+-f{k,;,ky)+ g/2], 

^{k^,ky,-K) = 2J^[l + -fika,, ky)] [1 - -fika,, ky) + g/2], 

where 

7(fca;, fcy) = ^(cosfcj, + COS/Cj^), (BIO) 


and g = J'/J. They satisfy the relation cuk = t^k+q, 
which implies that the dispersion of the symmetric and 
antisymmetric modes are identical up to a shift of the 
in-plane momentum by (7r,7r). The symmetric and anti¬ 
symmetric modes soften respectively at k = 0 and k = Q. 

We next calculate the single-mode contribution to the 
dynamical spin structure factor in terms of the Bo¬ 
golyubov quasiparticles. To this end, we consider the 
operator 




1 

^/Ns 

1 

7^ 


I 




+ ®k+Q 


-k-Q 



(Bll) 


In terms of the Bogolyubov quasiparticles, the effect of 
onto the vacuum |0) is given by 


‘S'k |0) = - Mk)aLk + (^^k+Q + Uk+Q)aL 


2 L 


k-Q 


| 0 ), 

(B12) 

where we have used the property Wk = w-k and Wk = 'C-k- 
Inserting Eq. (B12) into the Lehmann representation of 
the dynamical spin structure factor, and summing over 
the single-particle eigenstates of Eq. (B 81 gives 


Ss{uj, k) = Zs(k)(5(w - Wk), 
with the single-mode amplitude 
S 


(BIS) 


'2^s(k) = - (|uk - Wkp + |wk+Q + Uk+Qp) , (B14) 


where we have used the identity Wk = Wk+q. Integrat¬ 
ing Ss{uj, k) over uj yields the integrated spectral weight. 


^^(k), which within linear spin-wave theory in the single¬ 
particle sector is thus given by the single-mode amplitude 
Z 5 (k). The integrated spectral weights for the symmet¬ 
ric and antisymmetric channels are thus given by 


2s{kx,ky,0) 


^S^kx^ ky^ tt) 


4S^Jil-jikx,ky)) 

^{k^,ky,0) 


(B15) 


4S^Jil-jikx,ky)+g/2) ^ 

,ky ,7r) 


Thus we recover the limiting behaviour obtained numer¬ 
ically: the amplitude of the symmetric mode at k == 0 
vanishes linearly, while the amplitude of the antisymmet¬ 
ric mode at the antiferromagnetic Bragg peak position 
k = Q diverges like ~ l/|k — Q|, namely, 

JQL 

Zs{k,k,TT) Ki ^ fc^-TT, (B18) 

respectively. 


2. Higher-order spin-wave theory 


As discussed in the main text, the agreement with the 
quantum Monte Carlo data is improved if the spin-wave 
expansion is extended beyond the harmonic order. The 
first step is to truncate the Holstein-Primakoff expansion 
to next-to-leading order. This amounts to keeping, in 
the expression of the spin operators in the rotated frame 
{x',y',z'), terms of order such that 


sr 


Sf 

sf 


2 

~2r 


(o; -f aj) 

(fl; - aj) 


S-Tli . 


1 


4v^ 

1 


4iV^ 


(n^a, +a|n,), 
(n,a; -ajui), 


(B19) 


Given the collinear nature of the classical state about 
which we perform the expansion, the terms beyond the 
usual harmonic contribution will be of quartic order in 
bosonic operators. To this level of approximation, we are 
left with a Hamiltonian 


R = iL(o) -H -p ifC) o(i/s) , 


(B20) 


where and are respectively the classi¬ 

cal energy contribution ~ the harmonic fluctuation 
Hamiltonian ^ S, and the collection of all quartic inter¬ 
action terms ~ 1. Lastly, 0(1/S) denotes all the remain¬ 
ing terms in the expansion which are of order I/S'" with 
a > 1 and which we neglect. Th e expression for the sum 
jjio) _|_ fj{ 2 ) jg given in Eq. (B4), while is obtained 























14 



values in the harmonic ground state. This yields an ef¬ 
fective harmonic contribution ~ 1 which adds up to 
the quadratic terms ^ S. Given the main objective 
of this calculation, which is just to demonstrate that the 
qualitative tendency is correct, we have not attempted to 
go beyond this simple treatment and to make the calcula¬ 
tion of the expectation values self-consistent. This more 
sophisticated and significantly heavier approach, known 
as self-consistent spin-wave theory, includes higher-order 
corrections in an approximate way and would probably 
further improve the agreement, but it would anyway not 
allow to describe the quantum phase transition. 

In Fourier space the sum i?eff = + H^g 


has the same structure as Eq. (|B4|), with 
Hes = -S{S + 1) 


2J+^]Ns 


+ S 


k 


/ t \ f A f Bf \ 
rk-«-kj j 


flk 

t 

-k 


a 


+ C, 


where C is an additional constan t of order 0{1) and 
Af,Bf are the coefficients of Eq. ( |B5[ ) , evaluated at the 
effective couplings Jeff and J^g (i.elT^^® = ^k( Jeff, Jeff) 
and B^ = B\^{Jes, J'fj)), which are given by 

T Til A —n\ , . f A'— n 

JeS — J ( 1 H - ^ - ) , Jeff — J I 1 


S 


S 


(B21) 

where n, A and A' are the following two-body averages 
computed in the harmonic approximation 

n = (a|a;), 


FIG. 12. (Color online) Effective coupling ratio Oeff as a func¬ 
tion of gior S = 1/2. 






(.kx,ky,k^) 


FIG. 13. (Color online) Plots of the dynamical spin structure 
factor at different ratios g = J' / J inside the antiferromagnet- 
ically ordered region for the spin-half Heisenberg model on 
the square lattice bilayer. The panels compare the quantum 
Monte Carlo (QMC) results to both the linear and 2nd order 
spin-wave theory results as well as to further rescaled effective 
coupling ratios, which for g = 2 and g = 2.522 are obtained 
as g^g = 1.8244(; and g^g = 1.95g respectively. 


A — (a/Of+d) for d e {±x, ±y} , (B22) 


A' = {aiaigg) forde{±z}. 

Note that all other two-body averages vanish, (a|a;+d) ~ 

0, and (oj ) = 0. The next-to-leading-order term of the 
spin-wave expansion thus has the effect of renormalizing 
the coupling ratio g = J'/ J via 


Jeff — 9 


S + A' — n 
S + A — n 


(B23) 


Depending on the values of the bare couplings J 
and J', the multiplicative coefficient ranges from values 


smaller than one to values bigger than one, as shown 
in Fig. [T^ This figure suggests that for ratios 5 > 1, 
spin-wave interactions yield an effective ratio ^eff which 
is larger than the bare value. We conjecture that this 
renormalization will be even bigger if higher-order spin- 
wave interactions are taken into account. For example, 
we find that the quantum Monte Carlo results at 5 = 2 
are rather well htted by the spin-wave dispersion if the 
effective ratio, used is even larger than the one ob¬ 
tained from this improved spin-wave theory, as shown in 
Fig. This figure also shows that for g = 2.522, an 
almost similarly good fit can be obtained. 
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Appendix C: Perturbation theory in 1/g 


For vanishing intralayer coupling J = 0, i.e. in the 
limit l/g = 0, the ground state of the bilayer spin system 
is given by the direct product of singlets on the interlayer 
bonds (rungs). We denote this state by [S') = |s)i, 

where |s)i is the singlet on rung i, 

|s)i = (|t)i,l|i)i.2 - |i)i,l|t)i,2) • (Cl) 


To treat the intralayer coupling J as a weak perturbation 
in the small 1/g regime, we separate the total Hamilto¬ 
nian H of the bilayer model in the interlayer and in¬ 
tralayer parts i/_L and Hy, which scale respectively with 
J' and J. \S) is thus the ground state of H±. Treating 
i/|| as a weak perturbation to the ground state to 
first order in 1/g is found to be 

|G5) = I*5)-^ E [\toh\toh+d (C2) 

i,dG{x,y} 

l)i+d 1 ^—1)2 |^l)z-t-d] 5 

where the state \tn)i\tm)i+d is the direct product of the 
triplets \tn) and \tm) on rungs i and i-l-d times the direct 
product of singlets on all other rungs (which are thus 
implicit in this notation). To study the dynamical spin 
structure factor, we consider the following local interlayer 
bond operators which relate to the antisymmetric {kz = 


) and the symmetric {kz = 0) modes: 



(C3) 

= ^(Ei+E2 )- 

(C4) 


These operators act on the eigenstates of a rung as fol¬ 
lows: 


st^\h)^ = o 

S^^\to)i = 0 

= Is)* 


Sf\s). = 0 

Sflhh = 0 

Sf\to)^ = \ti)i 
Sf\t-l)^ = \to)i 


(C5) 


We first consider the antisymmetric sector. We start by 
computing the effect of 


1 

^/N 


^ \ ^ qAS 




(C6) 


on the approximate ground state Eq. (C2) (where k is 


a two-dimensional vector and N is the number of inter¬ 


layer bonds). According to (C5) this operator promotes 
a singlet to a triplet, and thus the antisymmetric mode 
captures the dynamics of triplet excitations. Restricting 
ourselves to the single triplon sector, we obtain: 

S^^\GS) = (^-1 + E(fc., ky)^ |ii)k , (C7) 


with 




(C8) 


and ^(kx, ky) defined in Eq. (BIO). To compute the ma¬ 


trix elements entering the expression of the dynamical 
spin structure factor in the Lehmann representation, the 
knowledge of the full spectrum is required. Here, we 
will consider the perturbative eigenstates. Since we con- 
side r on ly the single triplet component of S^^\GS) (see 
Eq. C7), we can restrict ourselves to the eigenstates of 


the effective Hamiltonian 


= H^+ Pii7||Pi -h PiH\\SHi \Pi (C9) 

where Pi is the projector onto the single triplet degener¬ 
ate manifold, and S = {1 —Pi)/{Eit —H±), Eu denoting 
the energy of the single triplet degenerate manifold. The 
states |ti)k are found to be eigenstates of 77®® with eigen¬ 
values 

{i-p'N + J') + 2J-f{kx,ky) 

-^(|(7V-2) + i(47(A:,„A:,)2-l)}|ti)k. 

(CIO) 

A similar calculation allows us to compute the ground 
state energy to the same order in J/J', 

Eo = -lfN-l^N. (Cll) 

Thus, overall we obtain for k = [kx, ky, tt)'^, i.e., in the 
antisymmetric channel, the single-triplon contribution to 
the dynamical spin structure factor as 


5 * 5 ( 0 ;, k) = Zs(k)i5(a; - Wk), (C12) 


with the spectral weight amplitude to first order in 1/g 
given by 

'2s(k) = 1 - ‘^j{kx,ky). (C13) 

9 

The single triplon dispersion, correct up to quadratic or¬ 
der, is given by 


Wk = J' + 2J^{kx,ky) - 2— ['^{kx,ky)'^ - l] . (C14) 


Equation ( C13| ) shows that the amplitude of the dynami¬ 
cal spin structure factor for the antisymmetric mode is of 
order 1. The momentum modulations of the amplitude 
are an effect of order 1/g, with the largest amplitude for 
k = (tt,TT,tt)''' and the lowest at k = (0,0 ,tt)'''. These re¬ 
sults compare well with the quantum Monte Carlo struc¬ 
ture factor, as discussed in the main text (see Fig. [^. 

Next, we turn to the calculation of the symmetric chan¬ 
nel. We start by computing, for k = {kx, ky^, the effect 
of 


1 

^/N 


^ \ ^ qS 




(C15) 









16 


2 J' + 4,7- 



FIG. 14. Spectrum of within the two-triplet sector for a 
system size L = 10. Energies are measured with respect to 
the ground state energy. 


on the approximate ground state in Eq. (C2). According 
to Eqs. (C5), this operator cancels the singlet component 
of the wave function. The finite contributions come from 
the triplet component, and are of order 1/g, 


^k|G5)=-f E (Mo)k.d - |Wi)k,d) (1 - , 

^ de{x,y} 

(C16) 

where we have introduced 

|^lio)k,d = E 6**^ *'‘l^l)*l^o)i+d , (C17) 

|^oil)k,d = E *'‘l^o)j|^l)i-l-d , (C 18 ) 

the Fourier transforms of neighboring pairs of triplets 
with zero relative momenta. Similarly as before, to com¬ 
pute the matrix elements required for the dynamical spin 
structure factor, we will consider the perturbative eigen¬ 
states. Since, to order 1/g, the states S'q|GS') contain 
two triplets, |to) and |ti), we can restrict the calculation 
to the eigenbasis of the effective Hamiltonian 

Hf = + (C19) 


value of the momentum, we computed the overlap be¬ 
tween S’^jGS') and the eigenbasis of Thus, for each 
momentum, there is a distribution of amplitudes at dif¬ 
ferent energies. This distribution of amplitudes is shown 
in Fig. [^of the main text. Based on our analysis, a few 
remarks can be made about the dynamical spin structure 
factor associated to the symmetric mode: (i) It can be 
seen analytically that at zero momentum the amplitude 
of the dynamical spin structure factor in the symmetric 
channel va nishe s (this may also be seen in Fig.™. In fact, 
from Eq. (C16) one obtains <S'^^g|G5') = 0. (n) At mo¬ 
mentum k = (tt, tt)''', one can verify that <S'k=(,r ,r)T 
is an eigenstate of with energy equal to 2J' — J12 
when counted from the energy of the groundstate. This 
explains why at A: = (tTjTt)''’ all the amplitude is concen¬ 
trated at a; = 2J' — J/2 in Fig. (hi) The amplitude 
of the symmetric mode is of order (1/g)^, while that of 
the antisymmetric mode is of order 1. Thus, the signal of 
the symmetric mode vanishes in the limit l/g —> 0 (see 
Fig.[l]). 


Appendix D: Bond-bond correlations from the 1/d 
expansion 

In this appendix, we provide details of the calcula¬ 
tion of the dispersion and the single-particle weight in 
the bond-bond correlations, which are related to the 
Higgs-peak contribution of the dynamical singlet struc¬ 
ture fact or 5'a (a;, k), using the bond-operator-based 1/d 

expansiorPSE], 

The starting point is a model of coupled dimers de¬ 
scribed by Eq. Q, but now on a hypercubic lattice in 
d spatial dimensions, such that the square-lattice bilayer 
corresponds to d = 2. To obtain a non-trivial large-d 
limit, it is crucial to properly scale the ratio between the 
interdimer coupling J’ and the intradimer coupling J. To 
this end we introduce the ratio 


where P 2 is the projector onto the two triplet degener¬ 
ate subspace. The effective Hamiltonian describes the 
two-body dynamics of two triplets (|to) and |ti)) in a 
sea of singlets. The full diagonalisation of idy inside the 
two-triplet sector is non trivial, and we performed it nu¬ 
merically. We considered a L = 10 lattice, yielding a two- 
triplet Hilbert space of dimension 9900. We performed a 
full diagonalisation, resolving the energies in momentum 
space. The spectrum of forms a continuum in the 
thermodynamic limit; for our lattice size we have access 
to 100 values of the momentum in the whole Brillouin 
zone, and for each momentum, we get 99 eigenstates with 
different energies. The spectrum, with energies relative 
to the energy of the ground state, disperses around 2J', 
the energy cost of promoting two singlets into triplets. 
It is shown in Fig. (14|. We do not obtain an explicit 
expression for the dynamical spin structure factor in the 
symmetric channel since the eigenbasis and eigenenergies 
of are obtained only numerically. Instead, for each 


q = Jd/J' (Dl) 

as our tuning parameter, with the notion that q is kept 
fixed upon taking the d —> 00 limit. We refer the reader 
to Ref. Iinifor an extended discussion of the 1 /d expansion 
approach and the general formalism. In the following, 
we will restrict ourselves to T = 0 calculations. Further, 
most explicit expressions will be restricted to the lead¬ 
ing order in 1/d. Note that to this order (i.e. at the 
harmonic level) of approximation, the quantum critical 
point is located at q = q^ = 1/2. 

As shown in detail in Refs. sni and sn the bond- 
operator approach yields a triply degenerate triplon 
(spin-1 excitation) mode in the quantum disordered 
phase (i.e. for q < q^) while in the antiferromagnetically 
ordered phase (i.e. for q > qc) one obtains two degenerate 
transverse modes (the Goldstone modes) and a longitu¬ 
dinal amplitude mode. In view of its further usage we 
quote here the leading-order result!^ for the dispersion 
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of the longitudinal mode (valid for g > 1/2): 




where k = (fci, /c2, kdY and 


7k = ^ ^ cos ki. 


(D2) 


(D3) 


For the interlayer bond-bond correlations we consider 
the susceptibility 

/ OO 

dte*‘"‘(TtSk(t)B_k(0)), (D4) 

-OO 

where Tt is the time-ordering operator, and 

(D5) 
(D6) 


Bi — Sii • Si2 , 
1 


Bk = 


^/N ■ 


By the fluctuation-dissipation theorem, the imaginary 
part of the susceptibility xb is proportional to the struc¬ 
ture factor Sb ds considered in the main text. Expressed 
in terms of bond operators (we follow the notation of 
Ref. SOI)) we obtain 


Bi — Aa^ia . ■ 


(D7) 


In the disordered phase, the tic (tj^) are annihilation 
(creation) operators for local spin-1 excitations, intro¬ 
duced within the bond-operator formulation. In the fol¬ 
lowing, we work in the single-mode approximation (see 
Ref. im for details), i.e., we will ignore multi-mode con¬ 
tributions to the bond-bond correlation function which 
would cause continua in the spectrum. As a result, we 
neglect terms containing more than two triplon operators 
in the operator product B\^B_\^. We are then lead to the 
following expression for the interlayer bond susceptibility 
in the disordered phase: 


xr(w,k)=fV 




|Jk,o<5(w), (D8) 


where i?2 = {t\ofia) up to 0(\/d). We thus do not obtain 
any spectral weight in this sector apart from a Bragg-like 
peak at a; = 0. (We note that a two-triplon continuum 
appears at order 1/d.) 

In the ordered phase, we performed the calculations to 
leading order in l/d, i.e., at the harmonic level. Here, we 
use generalized triplon operators i and P, obtained from 
the t and P operators after a suitable rotation in Hilbert 
spac^^. Note that the p and iy operators refer to the 
two degenerate transverse modes, while the tz operator 
corresponds to the longitudinal mode. Furthermore, we 
define the condensate parameter 


corresponding to the condensation of one of the triplon 
modes in the ordered phase. We then obtain 


Bi — Sii • S2i 


O ~..i. — 

—__|_ /' f • _|_ /' f 


iy W 


1 


(tjj, + iiz) 


1-f A2 I 

I-f A2 


(DIO) 


where Q = (tt, tt, ...)''' and in the last equation we have 
used the singlet projector Pi = 1 — (with a = 

x,y,z). This gives 


b^ = ^Pn 


A2 


4 l + X^ 


1 

iq,ziq-k.z + 


<5k,0 + ^2 (^k-Q,z + ^-k-l-Q.z) 

(Dll) 

^q,a;^q-k,a; + Pq^ytq-k.y 


1-f A2 


As noted above, we now ignore the last term in the above 
equation since it will either contribute to the contin¬ 
uum or give a 1/d correction to the single-particle peak. 
Hence, we do not obtain a single-mode contribution from 
the transverse modes. We are thus left with 


X^^"(a;,k + Q)=iV 


A2 


■ 3 A2 
4 1 + A2 


(5k.o^(t^) (D12) 


+ ^ _l_ ^2 [SN,zi^i k) -|- GN,zi~ai, k) 
+Ga,z{<^, k) -|- Ga,z{—^, k)] 

.\2 1 ^ 

= N -- + Z -rzT (5k,0<5(w) 


4 l-kA2 


(D13) 


A2 


+ Y7p^(uk,z + Vk,z)^ 


1 


1 


OJ - Wk.z W + Wk,2 


where we have used the fact that k resides within the 
first Brillouin zone, and the fact that k ± 2Q is equiva¬ 
lent to k. In the above expression, Gn and Ga denote the 
zero-temperature normal and anomalous Green’s func¬ 
tion, respectively. Also, Uk,z and Uk,z are Bogoliubov co¬ 
efficients, used in diagonalising the Hamiltonian to lead¬ 
ing order in 1/d, and at this point we refer to Ref. [41] for 
their explicit expressions. We can now read-off the mode 
weight corresponding to ajk,z from the above expression, 


Znik + Cl) 


A2 

1 + A2 


(Uk.z + Vk,z)^ 


2q-l 

2\/4g2 + ’ 

(D14) 


A = 


2g- 1 
2g + l ’ 


(D9) 


As discussed in Sec. HI, we need to account for the 
shift in the location of the quantum critical point within 
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(fc,*,) (t.tj) 


FIG. 15. Spectral weight distribution in the dynamic singlet structure factor for different values of g, as obtained from the 
bond-operator theory using the mapping in Eq. 


the bond-operator approach from the actual value, when 
comparing the above results to the quantum Monte Carlo 
data. In Fig. we show the single-mode contribu¬ 
tion to the spectral weight atop the corresponding dis¬ 
persion relation for different values of g, using the map¬ 
ping in Eq. Q to the corresponding values of q used in 
the above formula. We recall that, in addition to this 
single-mode contribution, a two-magnon continuum ap¬ 
pears at leading order in 1/d which we have not deter¬ 
mined here. Finally, we note that within the current 


approach, one can similarly calculate the quasi-particle 
contribution to the bond-bond correlations among the 
intralayer bonds. Similar to the case of the spin-spin 
correlations, one can in this case distinguish between a 
symmetric and antisymmetric channel with respect to 
layer inversion symmetry. For the symmetric channel, 
we obtain again a single-mode contribution like in the 
interlayer case considered explicitly above, however of 
reduced spectral weight. In the antisymmetric channel, 
no single-particle contribution is obtained so that only a 
continuum contribution results. 
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